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Abstract. Using a fundamental measure density functional theory we investigate 

both bulk and inhomogeneous systems of the binary non-additive hard sphere model. 

rH ■ For sufficiently large (positive) non-additivity the mixture phase separates into two 

O ! fluid phases with different compositions. We calculate bulk fluid-fluid coexistence 

curves for a range of size ratios and non-additivity parameters and find that they 
compare well to simulation results from the literature. Using the Ornstein-Zernike 
1,.^ _ equation, we investigate the asymptotic, r — > cxd, decay of the partial pair correlation 

ly-v ' functions, gij{r)- At low densities there occurs a structural crossover in the asymptotic 

^^ , decay between two different damped oscillatory modes with different wavelengths 

^^ ' corresponding to the two intra-species hard core diameters. On approaching the 

fluid-fluid critical point there is Fisher- Widom crossover from exponentially damped 
oscillatory to monotonic asymptotic decay. Using the density functional we calculate 



^— s ' the density profiles for the planar free fluid-fluid interface between coexisting fluid 

phases. We show that the type of asymptotic decay of gij{r) not only determines the 
asymptotic decay of the interface profiles, but is also relevant for intermediate and 
even short-ranged behaviour. We also determine the surface tension of the free fluid 
interface, finding that it increases with non-addivity, and that on approaching the 



C^ ' critical point mean-field scaling holds. 

PACS numbers: 61.20. Gy 64.75. Gh 68.05.-n 
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1. Introduction 

Liquids can consist of a mixture of components which in molecular systems may be 
different atomic or molecular components. In colloidal systems mixtures may be 
formed from particles with differing shapes or sizes, or mixtures of colloids and non- 
adsorbing polymer. Besides gas-liquid phase separation such mixtures may exhibit 
liquid-liquid separation, where the system demixes into two (or more) phases with 
differing compositions. 

Inter-particle interaction potentials between the constituent particles in a gas or a 
liquid may contain both short-range repulsion, due to the overlap of outer electron shells, 
and longer-ranged attractive or repulsive tails, due to dispersion or Coloumb forces [1]. 
For intermediate densities both features of the potential are important, and as van der 
Waals discovered p], it is the presence of an attractive tail that drives liquid-gas phase 
separation. However, if the liquid is dense, the long-range tail becomes less important 
and the structure of the fluid is primarily determined by the short-range repulsion. 

Following van der Waals [2] it is convenient to separate the potential into its 
short-range sharply repulsive and longer-range components and treat them within a 
theoretical approach separately. The simplest model for the short-ranged repulsive 
potentials is the hard sphere model which disallows particle overlap. This model 
has been shown to give a good approximation to the thermodynamic and structural 
properties of fluids, particularly near crystallisation. The development of ever better 
approximate theoretical treatments of the hard sphere model is a major element of liquid- 
state theories. Furthermore, given a theoretical treatment of the hard sphere potential, 
attractive or repulsive tails may be incorporated using relatively simple perturbation 
theories. It is straightforward to generalise the hard sphere model to multi-component 
mixtures. This simplest mixture has two components where the second species can have 
the same or a different diameter than the first species. This two-component mixture 
is normally formulated such that the distance of closest approach between particles of 
different species is a simple mean of the diameters of the particles of each species. In 
a real fluid this assumption may not be true and relaxing this constraint allows novel 
features of real systems to be investigated [SI SI IS] • 

Specifically, we define a two-component mixture of particles that interact through 
the hard sphere pair-potentials, 

{ oo r < ffij 
(1) 
otherwise, 

where i, j = 1, 2 label the species, and an are the particle diameters. The usual additive 
cross-species range of interaction is 0"i2 = |(o"ii + 0-22) but the non-additive hard sphere 
(NAHS) model generalises this so that ai2 can be smaller or larger than the arithmetic 
mean of the like-species diameters, 

c^i2 = -(l + A)(aii + a22), (2) 
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where A > — 1 measures the degree of non-additivity. We characterise the model by 
the size ratio, q = crii/cr22 < 1? and by A. For A = the model reverts to the binary 
additive hard sphere model. 

If both A and the density of the fluid are sufficiently large, the fluid demixes 
into two phases, one rich in particles of species 1, one rich in particles of species 
2 O [HI O |8l |9] . Experimental work on a number of systems, including alloys, aqueous 
electrolyte solutions, and molten salts, suggests that non-additivity may lead to both 
hetero-coordination and homo-coordination |5l|10]. In recent work Kalcher et al have 
used Monte Carlo simulations to calculate an effective interaction between charged ions 
in a electrolyte solution [TT], integrating out the degrees of freedom of the solvent 
molecules. Using a Barker-Henderson mapping they identified hard core diameters with 
values of non-additivity as large as 0.36 for Nal. In contrast, for A < mixing of the 
two species is encouraged and the fiuid can exhibit strong short-range order [12]. For 
A being sufficiently negative, clustering effects, mesoscopic ordering and the formation 
of heterogenous structures were reported [T2|, [131 [El- 

Besides the additive model there are two further important limiting cases of the 
binary NAHS model. The Asaka-Oosawa-Vrij model [151 [13 [E] has long been used 
as a simple description for the behaviour of a mixture of colloids and non-adsorbing 
polymers. The colloids interact through hard sphere interaction with diameter cTcc, 
while the polymers are treated as ideal, cTpp = 0. The colloid-polymer interaction is also 
hard sphere like but with a diameter acp = CTccf^ + Rg, where i?^ > is the polymer 
radius of gyration. In the formulation of ([2|) this corresponds to A = 2Rg/acc- The 
second important case is the binary Widom-Rowlinson model [IB] where both same- 
species interactions are ideal, an = 0-22 = 0, but the cross-species interaction is through 
a hard-core with diameter a^ > 0. This corresponds to taking the limit of ([2]) where 
A —7- 00 and (Xu = o"22 — )■ while keeping a^ constant. The Widom-Rowlinson model 
has become an important model in statistical physics due to its entropy driven demixing 
transition, yet simple structure of interactions. For a more thorough review of the 
literature on NAHS mixtures we refer the reader to the review article [5] and to the 
more recent contributions [T9l [20l[2Tl [22]. 

In this paper we explore the binary NAHS using the recently developed fundamental 
measure density functional theory for this model [23]. Fundamental measure theories 
(FMTs) are a class of density functional theories (DFTs) that are based on the 
geometrical quantities (fundamental measures) of the particles involved, e.g. volume, 
surface area, radius, and Euler characteristic. These fundamental measures enter the 
theory through weight functions which are based on these measures. The weight 
functions are convolved with the density profiles in order to give a set of weighted 
densities which are then combined within a free energy density. The original functional 
was formulated by Rosenfeld for additive mixtures of hard spheres [24j using both 
scalar and vectorial weight functions. Subsequently Kierlik and Rosinberg constructed a 
functional for hard spheres using only scalar weight functions [25]. Their theory was later 
shown to be equivalent to Rosenfeld's formulation [20] . These functionals reproduce the 
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Percus-Yevick approximation for the two-body correlation functions in the bulk fluid. 
One drawback of these original hard sphere FMTs was that they were not suitable for 
studying crystallisation phenomena due to unphysical divergences within the functional 
when the density profiles became strongly confined. This was later remedied (for hard 
spheres), first using a simple modification |2^, and then later by introducing tensorial 
weight functions p8]. For recent reviews on DFT and hard body DFTs in particular, 
we refer the reader to [221 EO]- In studies that preceded the NAHS functional, FMTs 
were proposed for both the Asakura-Oosawa [Blj and Widom-Rowlinson |32] models. 
There is a number of publications dedicated to the study of interfacial properties of the 
Asakura-Oosawa model |33l |3ll |35] . Note that the present DFT reduces to that used 
in [311 ES], when the Asakura-Oosawa limit of the general non-additive hard sphere 
mixture is taken. 

The NAHS functional, which was first introduced in [23] is based on the scalar 
Kierlik-Rosinberg deconvolution [25] of the hard sphere weight functions, but introduces 
a further ten scalar weight functions to take account of the non-additivity. It was shown 
that the functional correctly predicts fluid-fluid phase separation and that the theory 
provides a reasonable prediction for the location of the critical point compared to existing 
simulation results. Using the Ornstein-Zernike (OZ) equation, i.e. by inverse Fourier 
transforming the analytic (Fourier space) total correlation functions, it was shown that 
the theory provides good account of the radial distribution functions, gij{r), as compared 
to Monte Carlo simulation results, though these are not the same as those obtained by 
the Percus-Yevick (PY) approximation and they violate the core condition gij{r) = 
for r < (Tjj, where gij{r) is the partial pair correlation function between species i and j. 
The non-additive functional has also been formulated for the one-dimensional version 
of the model, binary non-additive rods on a line [36], making accurate predictions for 
the particle correlation functions, although failing to reproduce the exact solution [57] . 
In further work [38j, the spherical and one-dimensional convolution transforms in the 
theory were investigated and shown to form an Abelian group. 

In more recent work, Ayadim and Amokrane [39] have used the functional of [23] 
to calculate the radial distribution functions via the Percus test particle route [IQ]. This 
involves introducing an external potential that represents a single particle fixed at the 
origin and numerically solving for the density profiles around it. The core conditions are 
automatically satisfied. These results for the radial distribution functions can then be 
compared to those obtained from the simpler OZ route, in order to assess the internal 
consistency of the functional. The authors of [51] found that gij{r) calculated via the test 
particle route exhibit small but clearly noticeable unphysical jumps that are not present 
in the results from the OZ route. They argue that these are not numerical artifacts, 
but are due to shortcomings in the construction of the functional, and suggest that 
the functional requires changes at a fundamental level to eliminate the occurrence of 
discontinuities. In the present investigation in planar (rather than spherical) geometry, 
we do not find unphysical kinks in density profiles. Both planar fiuid-fiuid interfaces, as 
well as the density profiles near a planar hard wall |1T] are free of artifacts. Furthermore, 
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the current study is dedicated to the intermediate and long-ranged behaviour of the bulk 
fluid pair correlation functions. We expect this to be largely unaffected by the jumps, 
which were found to occur at short separation distances ^39| . 

For fluids where the pair-potentials are short-ranged, it can be shown that hij{r) = 
Qijir) — 1 can be evaluated by determining the positions and residues of the poles 
(divergences) of the complex structure factors, Sij{k), where fc is a complex wave- 
number [12]. These poles either occur as a complex conjugate pair, which give rise to a 
damped oscillatory contribution to rhij{r), or as a single purely imaginary pole, which 
gives rise to a purely exponentially decaying contribution. In general, there is an inflnite 
number of poles. However, in order to flnd the intermediate and asymptotic, r — )■ oo, 
decay of hij{r), it is usually sufficient to flnd the positions of a small number of poles 
- those that give rise to the slowest decaying contributions. Furthermore, the ultimate 
asymptotic decay is determined by the pole(s) with the smallest imaginary component, 
referred to as the leading order pole(s). As the model parameters (or statepoint) are 
varied the identity of the pole with the smallest imaginary component may change, 
leading to abrupt changes in the type of asymptotic decay. Furthermore, there may 
also be crossover in oscillatory asymptotic decay with different wavelengths. A number 
of studies have shown that changes in the asymptotic decay mode may be detected in 
simulation and experiments and have verifled its relevancy in studying the microscopic 
properties of the fluid [H Si SSI HHl HZ] . 

The raison d'etre of DFT lies in its prowess to investigate inhomogeneous situations 
in equilibrium. Given an approximation for the density functional, taking the derivative 
with respect to the density proflles yields a set of Euler-Lagrange equations. By 
numerically solving these equations one obtains the set of equilibrium density proflles, 
which minimise the grand potential functional of the system. 

In the present paper we consider the NAHS model with A > and calculate fluid- 
fluid demixing binodals and spinodals for a range of size ratios, g, and non-additivity 
parameters, A. For size ratio g = 0.1, we compare the predictions for the coexistence 
curves to results of computer simulations by Dijkstra [7]. We flnd that the theory 
reproduces the location of the binodal reasonably well. By calculating the partial 
pair direct correlation functions from functional derivatives of the excess free energy 
functional and inverting the OZ equation in Fourier space, we determine the positions 
of the poles, as well as the asymptotic, r — )■ oo, decay of the partial pair correlation 
functions. We flnd that besides structural crossover between oscillatory decay with one 
wavelength to oscillatory decay with a different wavelength (that also occurs in the 
additive model ^8]), there is also Fisher- Widom (FW) crossover from exponentially 
damped oscillatory decay to monotonic exponential decay. We use the functional for 
investigating the planar fluid-fluid interface between coexisting phases and demonstrate 
how the different types of asymptotic decay of the bulk correlations in the two coexisting 
phases determine the asymptotic and intermediate decay of the density proflles on both 
sides of the fluid- fluid interface, consistent with the general theory of asymptotic decay 
of correlations [12]. The essential quantity in studying fluid interfaces is the surface 
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tension, 7, which is the excess free energy per unit area required to maintain the surface. 
This can be measured experimentally and therefore provides a direct connection between 
theoretical approaches and real fluids. We determine 7 quantitatively and show that as 
the critical point is approached, mean-field scaling is reproduced. 

The paper is structured as follows: In section |2|, which can be safely skipped by 
expert readers, we outline the relevant theory including a description of the excess 
free energy functional used, the Ornstein-Zernike equation for binary mixtures and 
the theory of asymptotic decay of correlations. We also describe how to calculate 
inhomogeneous density profiles. The main results of the work are described in section [21 
where we present bulk phase diagrams for a range of size ratios and non-additivities. We 
investigate the pole structure of the structure factors in the complex plane and indicate 
the regions of the phase diagram with different types of asymptotic decay. We calculate 
the density profiles for the fluid-fluid interface and the surface tension. In section H] we 
present a discussion of the results. Finally, in the appendices we present further details 
on the structure of the weight and kernel functions, and explicitly write down the free 
energy density contributions. 

2. Theoretical Background 

2.1. Density Functional Theory 

We first reintroduce the general density functional theory framework |l9l [50]. For a 
classical system composed of two different species of particles, one can construct a grand 
potential functional VL[pi, p2\ of the set of one-body density profiles Pi(r) for i = 1, 2, 

2 r 
fi[pi,P2] = F[p„p,]-Y, drp,(r)(p, - Vr\r)), (3) 

j=i -^ 
where F[pi,p2] is the intrinsic Helmholtz free energy functional, pi is the chemical 
potential of species i, V^'"'*(r) is an external potential that acts on particles of species i 
and r is the spatial coordinate. The intrinsic Helmholtz free energy functional may be 
separated into two contributions: 

F [pi , P2] = Fid [pi , P2] + Fex [pi , P2] . (4) 

The first term in (jl]) is the Helmholtz free energy of an ideal gas, 

2 r 

i^id[pi,P2] = E^^^ / drA(r)(ln(A3p,(r)) - 1), (5) 

where ksT is the thermal energy and Aj is the thermal de Broglie wavelength of particles 
of species i. The second term in (jl]), Fex[pi,p2], is the excess contribution to the free 
energy which is due to inter-particle interactions. This part is in general unknown and 
is specific to the form of the inter-particle interactions. The FMT approximation for 
the binary NAHS excess free energy functional will be defined below. 

It can be shown that when fi[pi,p2] is minimised w.r.t. the density distributions, 
its value is equal to the thermodynamic grand potential of the system, Q, and that the 
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set of density profiles that minimise VL[pi,p2\ are the set of equihbrium density profiles, 
Pj(r). One can summarise these two statements as 

5Vt[pi,p2] 



^Pii.^) 



= 0, fi[pi,p2] = fi, (6) 

Pl:P2 

where the left hand side of the first equation represents the functional derivative of 
^[Pi;P2] with respect to Pi(r), i = 1,2, evaluated with the equilibrium density profiles, 
pi(r) and p2(r). 

Properties of bulk fluid states can be obtained from evaluating (^ with constant 
bulk densities, Pj(r) = p^, such that the Helmholtz free energy in bulk is 

Hp\, p'2) = F[p\, p',]. (7) 

The chemical potentials and the pressure are given respectively by, 

^^^{pl,p2) = y Q-b , P{PUP2) = y + 2^PiPi, (8) 

where V is the system volume. For a system that exhibits phase separation, coexistence 
curves (binodals) are obtained by finding pairs of statepoints for which the chemical 
potentials and the pressure are the same in the two phases, labelled A and B, i.e., by 
solving simultaneously the three equations; 

pW=p(^) and p^^=pf\ 1 = 1,2, (9) 

where pi is the chemical potential of species i in phase A and P^^'^ is the pressure 
of phase A (and similarly for phase B). The limit of mechanical stability of the system 
(spinodal) can be obtained from the (numerical) solution of det {d'^{J^/V)/dpidpj) = 0. 
At the spinodal the compressibility of the fluid becomes infinite and the correlation 
length of the fluid, which is the length-scale on which the fluid correlations decay, 
diverges to infinity. The binodal and the spinodal meet at the critical point, where 
the correlation length in the coexisting phases becomes infinite and the order parameter 
that describes the difference between the two phases vanishes, i.e. the two phases become 
indistinguishable from each other. 

2.2. The Binary Non-Additive Hard Sphere Excess Free Energy Functional 

In FMT the excess Helmholtz free energy functional is constructed from a set of weighted 
densities, nj (r), which are formed by convolution of the bare density profiles, pj(r), with 
a set of geometrically inspired weight functions, Wu (r), appropriate for the model. The 
index z/ labels the type of weight function. In the binary NAHS functional the weight 
functions are spherically symmetric so that the weighted densities are given by 



n: 



«(x)= / drpi(rK(|x-r|,i?i), z = l,2, (10) 



where z/ = 0, 1, 2, 3, and Ri = aii/2 is the hard sphere radius of species i. The weight 
functions represent the 'fundamental measures' of the model in question, i.e. their 
volume (z/ = 3), surface area (2), integral mean curvature (1), and Euler characteristic 
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(0). They have dimension (length)'^"^, and therefore the weighted densities also have 
dimension (length) '^~^. We define components of a free energy density that depend on 
the weighted densities, 

<i..,(x,x') = $„, (KHx)},K^Hx')}) , (11) 

where a,/3 = 0,1,2,3. These terms have dimension (length)"'*'''"^ and their full form 
is given below. The excess (over ideal) Helmholtz free energy functional, Fex[pi,p2], is 
given by a double integral over space and double sum over the geometric indices, 

^^^^ =11 I /dxdx'<|.„,(x,x')K.,(|x-x'|), (12) 

where the convolution kernels, Kap[r), control the range of non-additivity between 
unlike components. The kernel functions are isotropic and are similar to the weight 
functions, although they depend on a new length scale, 

i?i2 = A(/?i + R2) = au - -(^11 + ^22), (13) 

which is the difference between the cross-species diameter and the mean particle 
diameter. (Note that in general Ru 7^ cru/^). The kernel functions Kai3 have dimension 
(length) ~"~'^, therefore the products ^ai3Kaj3 have the correct dimension (length)"^, as 
required by flT^ . 

We use the (fully scalar) Kierlik-Rosinberg form for Wv{r, R). Hence the four weight 
functions used in Eq. (ITU]) are defined as 

Wz{r,R) = sgn(i?)e(i?-r), 

W2{r,R) = 5{R-r), (14) 

on 

Wo{r,R) = -^6"{R-r) + ^6'{R-r), 
671 iTxr 

where r = |r|, R = i?i,i?2, -R12, sgn(-) is the sign function, 9(-) is the Heaviside step 

function, 6{-) is the Dirac distribution, and the prime denotes the derivative w.r.t. the 

argument. Although the hard sphere radii Ri are strictly greater than zero, the factor 

sgn(i?) is included in order that the weight functions may be reused below for the kernel 

functions, where R12 may be negative (if A < 0). 

The set of convolution kernels are symmetric w.r.t. exchange of indices, Kai3 = Kf^a, 

so that there are only ten independent weight functions. Four of these are given via (TT^ 

with R = R12; 

Koo{r) = W3{r,Ri2), Koi{r) = W2{r,Ri2), (15) 

Ko2{r) = wi{r, Ru), Kos^r) = Wo{r, R12). 

The set of further weight functions, suppressing for notational convenience the 
arguments of wi/ (r, -R12), are: 

Kii{r) = w[ = sgn{R)6'{R — r), 
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K22ir) 



W'n 



1 

Stt 
1 



5"{R 



w 



-1 



Ki3{r) = w-i 

K22,{r) = W_2 



647r2 
sgn{R) 
1 



(5(3) (i? 



—6"(R-r)-—6^^\R-r) 
27rr ^ ^ Stt ^ ^ 



(16) 



i^33(r) 



w. 



IGvr^r 
sgn(i?) 



(5(=^)(^-^)-7:TT'^^'n^ 



577^ 



647r2 
-6'-^\R-r) + l6^''\R- 



= Ri2, and the derivatives of the Dirac delta function are defined by 
(r'6{x)/dx'^ for 7 = 3,4,5. The Fourier space expressions of all weight 



where R 

<5W(a;) = 

functions are given explicitly in Appendix A 

The terms $q,^ are built from a sum of derivatives of the zero-dimensional excess 
free energy, 4>od{v) = {^ ~ v) 1^(1 ~ v) ~^ V^ where ?7 is a dummy argument (which can 
be viewed as the average occupation number of a zero dimensional cavity [27|), and 7 
labels the derivative: 0^ (?7) = (P (j)odiv) / '^V^ ■ The derivatives of 0od(^) are multiplied 
by products of weighted densities to ensure the correct dimensionality of the free energy 
density. We introduce ansatz functions Aa^ that possess the dimension of (length) ""^ 
and the order 7 in density (i.e. they contain 7 factors nj ). These are combined as 

*.. = EE«(U-,K?(4" + „f). 

7=0 7'=0 

Expressions for the non- vanishing terms of the ansatz functions are. 



(17) 



^01 



n, 



A 



(i) 



11 



n 



(i) 
' 

(i) 



|W 



W^« 



(0 



^02 ~ "-1 "■2 5 ^03 



A 



(i) 



"■1 "-2 : 

^Jn^2^] 



A 



(i) 



h^ 



2 J ; 



n 



(i) 



A 



(i) 



^12 ~ 8^V"'2 ; 5 ^21 ~ "-2 5 ^30 ~ -'-• 

The specific form of flTSl) ensures both that the terms in the sum in flT2|) possess the 
correct dimension of (length)"^ and that the prefactor of 0od in flTT|) is of the total order 
7 in densities. 

This completes the prescription for the excess Helmholtz free energy functional. 
Evaluating the sums in Eq. (I17p explicitly results in a total of 49 terms, which can 
be grouped either by the 16 kernel functions, or alternatively by the 10 unique weight 
functions. In Appendix B we transcribe some of the terms; all further terms can be 
obtained by symmetry. 



2.3. Fluid Structure and Asymptotic Decay of Correlations 

In order to study the pair structure of the bulk fluid, rely on the OZ equation, which 
separates the partial pair distribution functions, gij{r), into a 'direct' part between pairs 
of particles, and an 'indirect' part that comes from the interaction between all the other 
particles in the system: 



hij{r) 



Ctj[r^ 



1=1 



dr' hii{r')cij{\r -r'\), 



(19) 
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where hij{r) = gij{r) — 1 is the total correlation function and Cij{r) is the two-body 
direct correlation function between species i and j. The latter can be obtained from the 
excess free energy functional via functional differentiation, 



Q,(|r-r'|) = -(A;^T)-i 



6pi{T)6pj{T') 

By Fourier transforming one can re-write flTOl) as 



(20) 

/3l,P2=const 



kj{k) = Ci,{k) + Y,pThii{k)cij{k), (21) 

1=1 

where hij{k) is the (three-dimensional) Fourier transform of hij{r), 

hij{k) = — / dr r sin(fcr)/ijj(r), (22) 

^ Jo 

and similarly for Cij{k). It can be shown by rearranging the OZ equations that 

^.,,.).||, (23, 

where the common denominator is 

Dik) = [1 - piCii(A;)][l - p2C22{k)] - pip2Cu{k), (24) 

and the numerators in (1251) depend on the species indices: 

iVii(A;) = cnik) + p2[cL(/^) - Cnik)c22ik)], 

N22{k) = C22{k) + pi[c\^{k) - cn(A;)c22(fc)], (25) 

Ni2{k) = c,2{k). 

Using the definition of the inverse Fourier Transform, we obtain 

1 f°° 

^iji^) = TT^r dkksm{kr)hij{k), 



27r2r7o ^^ 



1 r°° N-(k'\ 

I dkksmikr)-^^^^-^. (26) 



2vrV Jo t){k) 

For the present functional, expressions for Cij{k) can be obtained analytically via (1201) 
and hence can be substituted into flM|) and fl25|) . before numerically Fourier transforming 
to obtain gij{r) = hij{r) + 1 from fl26|) . No numerical scheme for solving flTOl) is required. 
Indeed this method has already been successfully used in [23] and [36] to calculate the 
distribution functions and partial structure factors. 



S,,{k)=6,, + Jp\p%,{k), (27) 



where 6ij is the Kronecker delta. 



Another method that we will make extensive use of in the following is to investigate 
the singularities of hij{k) in the complex /c-plane [12]. Using ( 12B]1 and assuming that 
the singularities of hij{k) for the present systems are simple poles, we can proceed via 
Cauchy's residue theorem. Performing contour integration around a semicircle in the 
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upper half of the complex fc-plane, the total correlation functions can be written as a 
sum of contributions from the poles enclosed, 

rhij{r) = J2 ^n^ exp(iA;„r), (28) 

n 

where n labels the poles, fc„ satisfies D{kn) = 0, An is the amplitude associated with 
the pole at kn and i is the imaginary unit. The amplitude is related to the residue -Rn 

The poles are either purely imaginary, fc„ = iao, or occur as a complex pair, 
kn = ±«i + i«0 5 where both a^ and ai are real. In general there will be an infinite 
number of poles and contributions from many of those are required to account for the 
behaviour of hij{r) at small distances r. However, the ultimate, r — )■ oo, decay of all 
hij{r) is determined by the pole(s) that gives the slowest exponential decay, i.e., the 
pole(s) with the smallest imaginary part oq. These are referred to as the leading order 
pole (or poles in the case of a conjugate complex pair). 

If the leading order pole is purely imaginary, then all rhij{r) ultimately decay 
exponentially, rhij ~ Aijexp{—aor), as r — )■ oo, where Aij is an amplitude specific to 
each correlation function. On the other hand, if the leading order poles are a conjugate 
pair, then the sum of contributions from this pair of complex poles gives damped 
oscillatory ultimate decay, rhij{r) ~ 2Ajj exp(— oq^) cos(air — 6ij), with a common 
characteristic decay length a^^ and wavelength of oscillations 27r/Q:i. Aij and Oij denote 
the amplitude and the phase, respectively. 

As the model parameters and statepoint change, the positions of the poles in the 
complex plane vary. The pole(s) which have the smallest imaginary part, referred to 
as the leading order pole(s), can therefore be replaced by a different set of poles. This 
can lead to abrupt changes in the type of decay, either between damped oscillatory 
decay and monotonic decay, or between damped oscillatory decay with one wavelength 
to damped oscillatory with a different wavelength. 

For the one-component hard sphere fluid the decay is always damped oscillatory 
with a wavelength similar to the hard sphere diameter pE] • For the binary additive hard 
sphere mixture there is an abrupt crossover in the phase diagram from one wavelength 
similar to the diameter of species 1 to a different wavelength similar to the diameter 
of species 2 [IH]. The two different oscillatory wavelengths are each described by a 
complex conjugate pair of poles with real components which determine the oscillatory 
wavelength. This abrupt change occurs when these two pairs of poles have the same 
imaginary component. This marks a structural crossover line in the phase diagram. 

In general, there may also be a crossover between damped oscillatory and monotonic 
decay, particularly in systems which exhibit phase separation and where correlation 
functions obey Ornstein-Zernike (asymptotic exponential decay) behaviour close to the 
critical point. This crossover can occur via two mechanisms: In Fisher-Widom crossover 
a pair of leading order complex poles and a single imaginary pole change their positions 
in the complex plane as the statepoint is varied. As the critical point is approached the 
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leading order pole(s) change from the complex pair to the purely imaginary pole. The 
statepoints where this crossover occurs traces the FW line in the phase diagram. 

It has been shown, both via simulation and theory, that additive mixtures of 
hard spheres with small size-ratio, q < 0.2, can exhibit (meta-stable) fiuid-fiuid 
phase separation [51]. However, the PY approximation is unable to account for this 
phenomenon. Correspondingly, the asymptotic decay of the distribution functions is 
always oscillatory, for all size ratios, within PY theory [48J. In the same study, the 
authors also consider an effective one-component depletion potential, from which they 
are able to obtain the phase transition and to also find Fisher- Widom crossover form 
oscillatory to monotonic decay. Since the NAHS functional, taken in the additive limit 
A = 0, recovers the PY approximation for the bulk correlation functions, we do not see 
Fisher- Widom crossover in the additive model. 

In Kirkwood crossover two purely imaginary poles come together, coalesce and 
become a pair of complex poles. This mechanism often occurs in fluids that interact via 
soft, steeply-repulsive, pair potentials [521 153]. Indeed we do not find it in the present 
system. 

2.4- Inhomogeneous Systems 

In order to calculate inhomogeneous density profiles in the grand-canonical ensemble 
we consider the thermodynamic grand potential functional ([2]) and minimise Q[pi,p2] 
with respect to variations in the density profiles. This is equivalent to solving a pair of 
Euler-Lagrange equations, 

fi, = knTlogiA^pi) - kBTcf\r) + ^^^(r), z = 1,2, (29) 

where c\ (r) = —{kBT)^^6Fex/5pi{r) is the one-body direct correlation functional for 
species i = 1,2. In practice, the pair of equations (12^ must be solved simultaneously 
via an iterative numerical scheme. Explicit functional derivation of the excess part of 
the present functional, flT^ . yields 

= ^ Z'dxt.Wdx-rl) 

« — n ^ 



5J-e: 



^^^(^^ ,=0 



3 



(30) 



J2 /dx'^(x,x')ir„,(|x-x'|) 

which has the structure of two nested convolutions. For each value of 7, the partial 
derivatives of the excess free energy terms, 0a^^ = d^afs/dn^ , are first convolved with 
the kernel functions, Kai3{r), and then the sum of these is convolved with the single- 
particle weight function, w^{r). There are a total of 60 terms of the form of (I5U]) to be 
evaluated. 

In the present study we consider the planar fiuid-fiuid interface between co-existing 
fluid phases, where V^'^^*(r) = 0. The boundary conditions are chosen so that the density 
profiles decay to the coexisting bulk values far away from the interface. 
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Figure 1. Bulk fluid- fluid binodals (solid lines) and spinodals (dashed lines) for the 
binary non- additive hard sphere fluid with fixed size ratio q — owj 021 = 1 and varying 
non-additivity parameter A = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5 and 1 (from top to bottom), 
plotted as a function of the partial packing fractions, 771 and 772. For each system the 
binodal meets the spinodal at the bulk critical point (•). 




0.30 



Figure 2. Same as figure [11 but for fixed size ratio q = (J\\I<J22 = 0.5 and varying 
A = 0.05, 0.1, 0.2, 0.3, 0.4 and 1 (from top to bottom). 



3. Results 



3.1. Fluid Demixing Phase Diagram 



It has previously been shown [23] that for suitable parameters the present theory 
reproduces the phenomenon that the mixture separates into two different fluid phases [6] . 
Using ([8]) and solving for coexisting states we have calculated the coexistence curves 
(binodals), as well as the spinodals, for a range of model parameters. The statepoint 
is specified by the partial packing fractions, 77, = np'^afjG; recall that p^ is the bulk 
number density for species i. 

Figure [U shows the binodals and spinodals for the symmetric mixture, q = o\\jo22 = 
1, with A varying between 0.05 and 1, in the {rii,ri2) plane. Increasing A causes phase 
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0.020 



Figure 3. Same as figure [TJ but for fixed size ratio q — 0.1 and varying A — 0.1, 0.2, 
0.3, 0.4, 0.5 and 1 (from top to bottom). Note the scale on the (horizontal) 771-axis. 



separation at increasingly lower densities, and therefore the partial packing fractions at 
the critical point, ?7f"*, both decrease monotonically as A increases. Figure |2] displays the 
binodals and spinodals for asymmetric mixtures with fixed q = 0.5 and A again varying 
between 0.05 and 1. On increasing A, again both rj"^^ and 773"* decrease monotonically. 

We also consider a mixture with large size asymmetry, q = 0.1, where we can 
compare to Gibbs ensemble Monte Carlo simulation results of Dijkstra [7J. This large 
asymmetry is a significant test for the functional, as previous studies have showed that 
FMT struggles with large size asymmetry already in the additive case jSH |S5]. Figure El 
shows the binodals and spinodals for g = 0.1 and A varying between 0.1 and 1. Although 
we show the binodal for A = 0.1, in simulations it was found that fiuid-fiuid phase 
separation for this value of A is metastable with respect to crystallisation [7] . Note that 
the relevant range of values of 771 is much smaller than in figures [1] and [2j 

To compare our results to those of [7], in figure S] we plot the binodals in the plane 
spanned by pressure, P, and relative concentration of the larger particles, X2 = P2/P, 
alongside the simulation results of Dijkstra and the results from Barboy and Gelbart's 
mean- field theory [SS] (data taken from [7]). We find that the coexistence curves from the 
three approaches have similar shapes and positions, for all values of A shown. However, 
both theories predict demixing pressures that are lower than the simulation results. 
Such a systematic error is often a feature of mean-field theories, which underestimate 
the strength of density fluctuations close to the critical point. 



3.2. Asymptotic Decay of Correlations 

Despite the presence of subtle artifacts [39] in the results for the radial distribution 
functions, Qijir), calculated via the test particle route from the present functional, it 
has previously been shown that there is good agreement between the gij{r), calculated 
via the Ornstein-Zernike route and Monte Carlo simulation data, both in 3D [23], and in 
ID [36]. Here we explore in detail the asymptotic, r — )■ 00, decay of correlations which 
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FM DFT critical points 
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Figure 4. Fluid-fluid coexistence curves for fixed size ratio q = 0.1 and varying A = 
0.2, 0.3, 0.4 and 0.5, plotted in the plane of pressure, P, and relative concentration of 
the large particles X2 = pij p- Results from present density functional (solid lines) are 
compared to the results from Gibbs ensemble simulations (crosses), as reported in [7], 
and a mean- field theory due to Barboy and Gelbart |56j (dashed lines). The pressure 
is scaled by the size of the larger species so that the figure is consistent with figure 3 
of [?]■ Both theories predict phase separation at lower pressures than the simulation 
results. 
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Figure 5. The partial radial distribution functions, gij{r), for the system with q — 1 
and A — 0.1. Note that gii{r) = g22{r) and that the pairs of curves are offset 
upwards by 4 units for clarity. The pairs of profiles are at the following statepoints: 
(?) r]i = i]2 = 0.1, {ii) i]i = ri2 = 0.13, and {Hi) r/i = ??2 = 0.14. As the density 
increases and the statepoint approaches the binodal, the oscillations in all gij (r) become 
more pronounced. The inset shows ln\r hij{7')\ where hij{r) = gij{r) — 1 is the total 
correlation function. For (i) the intermediate and asymptotic decay is oscillatory with a 
wavelength ^ an. For {ii) the decay is oscillatory at small r, but as distance increases, 
the relative amplitude of the oscillations quickly decreases, and the profiles start to 
decay monotonically. For {Hi) there is some short-range oscillations which die out by 
r ~ Tail so that the intermediate (and asymptotic) decay is monotonic. 



is determined by the poles of the partial structure factors, Sij{k), in the complex plane. 
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Figure 6. Same as figure [SI but for q = 0.5 and A ~ 0.1. The sets of profiles 
correspond to the following statepoints (i) rji = 0.151, 772 = 0.001, {ii) 771 = 0.051, 
?72 = 0.051, and {Hi) -qi = 0.151, 772 = 0.07. In the three subplots we display In \r hij{r)\ 
for the three different statepoints. In (i) the intermediate and asymptotic decay is 
oscillatory with a wavelength ^ an. In {ii) the intermediate and asymptotic decay 
is oscillatory is with a wavelength ~ (722 = San. In {Hi) the intermediate decay is 
oscillatory, but the asymptotic decay is monotonic. 



We expect the asymptotic decay to be robust and not be affected by tlie test particle 
artifacts. 

We begin by showing representative examples of the radial distribution functions 
obtained from numerically Fourier transforming the analytical expressions for 
hij{k) (12^ . For the symmetric mixture, q = 1 and A > 0, by varying the statepoint 
and calculating gij{r), we find that there are two types of intermediate and asymptotic, 
r — )■ cxD, decay. In figure Owe plot gij{r) for the symmetric mixture, q = 1, with A = 0.1, 
at three statepoints (i) r/i = 772 = 0.1, {ii) r/i = 772 = 0.13, and {Hi) rji = ri2 = 0.14. 
We find that as the statepoint approaches the binodal, the oscillations in gij{r) become 
more pronounced, when viewed on a linear scale. 

To elucidate the intermediate and asymptotic, r — )■ 00, decay of gij{r), the inset of 
figure shows In \rhij{r)\. Recall that hij{r) = gij{r) — l is the total correlation function. 
For the low density case, {i), the intermediate and asymptotic decay is oscillatory with 
a wavelength ~ an. As the statepoint approaches the coexistence region, the decay 
of gij{r) starts to become monotonic. For {ii) the intermediate decay is oscillatory, 
but the amplitude of these oscillations quickly decreases with increasing distance r and 
the decay becomes monotonic. For {Hi) the oscillatory contribution decays much more 
rapidly and the asymptotic and even intermediate (r > Tern) decay is monotonic. 

Calculating gij{r) for the asymmetric case g < 1 with A > we find that there 
are three types of intermediate and asymptotic, r — )■ 00, decay. In figure |6] we 
plot representative examples of gij{r) corresponding to the three types of decay for 
parameters g = 0.5 and A = 0.1. In the main panel of figure Owe show the set of 
gij{r) for the three statepoints {i) rji = 0.151, 772 = 0.001, {ii) 771 = 0.051, 772 = 0.051, 
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Figure 7. Parts (a) and (b): Complex partial structure factors, Sij{k), where k is the 
complex wave-number, for the symmetric mixture, (7=1, with A = 0.1, at statepoint 
771 = ?72 = 0.12. The height and colour of the surfaces represents the absolute value 
and the complex argument, respectively, of Sij{k). The surfaces are plotted over the 
range A: = to fc = 20+ lOi. Since the mixture is symmetric the intra-species structure 
factors are identical, i.e ^n (fc) = 522 (fc)- Note that the two structure factors both have 
strong divergences (poles) at identical positions. Part (c): The real partial structure 
factors Sij{k) plotted along the real axis, at the same statepoint as (a) and (b). 



and (Hi) rji = 0.151, 772 = 0.07. On the linear scale there is relatively small variation 
in the overall magnitude of the correlation functions between the three statepoints. In 
the subplots of figure |6] we display In \rhij{r)\ for the same statepoints. For statepoints 
close to the r/i-axis (i), we find that the decay is oscillatory with a wavelength ~ an. 
For statepoints which are close to the ?72-axis (ii), the intermediate and ultimate decay 
is oscillatory with wavelength ~ cr22. Since cr22 = 2o"ii we find that the oscillatory 
wavelength in (ii) is approximately twice as large as that in {i). Again, as we approach 
the coexistence region {Hi), we find that the intermediate decay is oscillatory, but the 
relative amplitude of the oscillations decreases with increasing r and the ultimate decay 
is monotonic. 

In order to understand this behaviour, we next determine the pole structure. For 
the one-component hard sphere fluid (g = 1, A = 0), it has been established that for 
all statepoints there is an infinite number of complex poles, fc„ = ai + iag, but there 
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Figure 8. Positions of tlie complex poles, fc„ = ai + iaoi of the partial structure 
factors, Sij{k) for q = 1 and A = 0.1. Only the poles with ai > are shown, and ag 
and ai are scaled by the diameter of the particles, an = CT22- The poles are labelled 
with an (arbitrary) index, n = to 6. The parts labelled (a) to (d) correspond to the 
points marked in the phase diagram, figure 111 (a) 771=772= 0.08, (b) 771 = 7/2 = 0.1, 
(c) 771 = 772 = 0.12, and (d) rji = r]2 = 0.14. As r] = r/i + ri2 increases, the imaginary 
component, ao, of the purely imaginary, 71, = 0, pole decreases and the leading order 
pole (the one with the smallest ao) changes from the pair of complex, 77 = 1, poles 
(n) to the single imaginary, 77, = 0, pole (•). The crossover occurs at the statepoint 
771 = 7/2 = 0.117. 



Binodal 

Spinodal 

Fisher-Widom line 

Tie-lines 




Figure 9. Fluid-fluid demixing phase diagram for q = 1 and A = 0.1. The tie-lines are 
at pressures Paf^/ksT = 3,5,7,9,11 (from bottom to top). The Fisher-Widom line 
(dash-dotted line) separates regions of the phase diagram where the decay of correlation 
functions has a different type (monotonic or damped oscillatory) of asymptotic, r — > 00, 
decay. 



are no purely imaginary poles [H]. Therefore, the asymptotic decay of the distribution 
functions, which is determined by the pole(s) with the smallest imaginary part, ao, will 
always be damped oscillatory, rhij{r) ~ 2Ajj exp(— aor) cos(air — %) and since g = 1 
there is only one length-scale in the fluid, so the oscillatory wavelength, 2-7r/ai, is always 



an. 
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Figure 10. Same as figure [71 but for q — 0.5 and A = 0.1, at statepoint rji = 0.111, 
7]2 — 0.051. The surfaces are again plotted over the range fc = to fc = 20 + lOi. Since 
the mixture is asymmetric, the intra-species structure factors, which are plotted in (a) 
and (c), are no longer identical. The real partial structure factors are shown in (d). 



For the symmetric mixture, g = 1, but with non-additivity A > 0, we find an 
infinite number of complex poles, as in the additive case, but there is also a single 
purely imaginary pole. In order to illustrate how the poles appear in the complex 
structure factors, in figure [7] we plot Suik) = 5*22 (fc) and Si2{k) as a function of the 
complex wave-number k, for parameters g = 1, A = 0.1, and statepoint rji = rj2 = 0.12. 
The height and the colour of the surface plots represents the amplitude and the polar 
argument, respectively, of Sij{k). Although Sii{k) and 6*12 (fc) are very different, they 
both exhibit sharp divergences at identical positions in the complex plane. These 
are the common poles of the complex partial structure factors. They are located at 
solutions of the equation, D{k) = 0, where D{k) is the common denominator fl2^ of 
the complex structure factors. To determine the positions of the poles, we numerically 
solve D{kn) = for complex fc„ = cti + iag. The relationship between the complex 
structure factor(s) and their more commonly known real structure factor(s) is that the 
former evaluated along the real axis equals the latter, see figure [7](c). 

Figure [8] displays a sequence of positions of the poles in the complex plane as 
the statepoint approaches the coexistence region, as indicated in the phase diagram in 



Non-additive hard sphere mixtures 



20 



5 
4 
t?3 
i5"2 
1 

5 
4 

l="2 

1 




(a) 






♦ . 


1/1=0 


□ 


V 
^ 3 


4 


- 


1 


2 


■ 


. 




▲ 


. 


■ 






■ 



(b) 












V 


♦ ■ 




D 




" 


■ 






" 



(c) 






■ 


1 


D 


V 


♦ " 


■ 






■ 



(d) 






■ 


■ 




V 


♦ 






▲ 




1 


n 















2.5 5 7.5 10 



2.5 



7.5 



10 



Figure 11. Same as figure HI but for q — 0.5 and A = 0.1. Tiie parts labelled (a) to 
(d) correspond to the points marked in the phase diagram, figure [T^ (a) rji = 0.151, 
7]2 = 0.001, (b) 7/1 = 0.051, ?72 = 0.051, (c) r?i = 0.111, 772 = 0.051, and (d) r^ = 0.151, 
7]2 = 0.051. In part (a) the leading order poles are a complex pair, n — 2, (a). In (b) 
and (c) the leading order poles are a different complex pair, n = 1, (n), and in part 
(d) the leading order pole is the purely imaginary pole, n — Q (•). 



figure [9l We show the positions of six of the complex poles, with index n = 1 to 6 
(arbitrarily) labelhng the poles, along with the single imaginary pole, n = 0. For each 
conjugate complex pair of poles, k = ±ai + iao we show only the pole with real part 
«! > 0. In part (a), at statepoint rji = rj2 = 0.08, the complex pair, n = 1, of poles with 
real part auai ~ 2ti are the leading order poles and give rise to ultimate oscillatory 
decay with a wavelength In ja\ ~ an. As we increase the total density, the positions of 
all poles change, but in general their imaginary components, ao, decrease, see figure[8]^b) 
for rji = ri2 = 0.1. 

The decrease in the value of a^ proceeds much more rapidly for the purely imaginary 
pole, n = 0, than for the complex poles and for r] > 0.117 this pole possesses the smallest 
imaginary part and therefore becomes the leading order pole. Figure IHI^c) is at statepoint 
^71 = ''72 = 0.12, where the leading order pole is the (n = 0) purely imaginary pole. 
This determines the ultimate asymptotic decay of correlations to be purely exponential. 
Increasing 77 further results in the imaginary components of all poles decreasing. In part 
(d), at statepoint rji = ri2 = 0.14, which is very close to the bulk critical point, the 
purely imaginary pole, n = 0, is still the leading order pole, and is very close to the 
real axis. As the critical point (or in general the spinodal) is approached the purely 
imaginary pole approaches the real axis, which corresponds to the divergence of the 
correlation length, 1/ao- 

By varying rji and r]2 we determine the statepoints where the leading order pole(s) 
changes from the n = 1 complex pair of poles to the n = purely imaginary pole. This 
yields the FW crossover line. We find that the FW line lies between the spinodal and 
the axes, and that the two ends of the FW line approach the spinodal - see figure M 
which displays the FW line alongside the binodal and spinodal for the case q = 1 and 
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Figure 12. Same as figure[Hl but for size ratio q — 0.5 and A — 0.1. Tire tie-lines are 
shown at pressures Pafi/kBT = 1, 1.25, 1.5, 1.75 and 2 (from bottom to top). There 
are three regions of the phase diagram, each with its own type of asymptotic decay. 
The structural crossover line (short-dash dotted line) separates the two regions with 
oscillatory decay and the Fisher- Widom line (long-dash dotted line) separates the two 
regions where the asymptotic decay is oscillatory from the region where the decay is 
monotonic. The structural crossover line for the additive case, A = (double dashed 
line), is shown for comparison. 



A = 0.1. Note that the FW hne intersects the binodal twice. Since the mixture is 
symmetric (g = 1) this occurs at coexisting statepoints. We show the importance of 
this feature below, when we investigate the planar fluid-fluid profiles. 

Turning to the asymmetric non-additive mixture, we find that there is again an 
infinite number of complex poles with non-vanishing real parts, as well as a single 
purely imaginary pole, see figure [TOl In figure [11] we plot the positions of the poles, 
with (arbitrary) index ra = to 4, at four different statepoints for parameters q = 0.5, 
and A = 0.1, as indicated in the phase diagram in figure [12j Figure [TTT a) plots the 
positions of the poles for statepoint rji = 0.151, 772 = 0.001, which is very close to the 
r^i-axis. The leading order poles are a complex pair, n = 2, that give rise to ultimate 
oscillatory decay with a wavelength, A = 27r/Q;i ~ 1.2crii. Figure ITlT b) is at statepoint 
?7i = 0.051, ri2 = 0.051 where the leading order poles are a different complex pair, n = 1, 
that gives rise to an oscillatory wavelength, A = 2.83o"ii. As the coexistence region is 
approached, the imaginary components of all the poles decreases, see figure(c), which 
is at statepoint rji = 0.111,7/2 = 0.051. This decrease in the imaginary components 
proceeds most rapidly for the purely imaginary, n = 0, pole which in figure ITlT d) 
{rji = 0.151, r72 = 0.051), now possesses the smallest value of ao and thus becomes the 
leading order pole. 

Therefore, for g < 1 and A > by varying the statepoint we find three regions 
of the phase diagram that have different types of asymptotic decay. There are two 
regions, close to the axes, where the decay is damped oscillatory with a wavelength 
similar to the majority component, and there is one region, which contains the spinodal 
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Figure 13. The positions, fc„ — ai+ iao, of the poles of the structure factors, Sij{k), 
in the complex plane for q = 0.5, r/i = 772 = 0.1 and increasing A from to 0.164. 
The positions of the poles are indicated for A = 0, 0.054, 0.108, 0.162. At A = 
there exists an infinite set of complex poles (no purely imaginary poles), four of which 
are shown here, n = 1,2,4 and 5 (+). The asymptotic decay is necessarily damped 
oscillatory, determined by the n = 1 pair of conjugate complex poles. As A is increased 
from zero, a second set of poles, three of which are shown n — 0, 3, and 6 (•), including 
one purely imaginary pole {n = 0), appears. Initially these poles have large imaginary 
components i.e. large ap. As A is increased, the imaginary components, ao, of the new 
set of poles decreases (the poles move down the complex plane). At A = 0.101 the 
leading order pole changes from the n = 1 pair of poles in the original set to the purely 
imaginary, n = 0, pole in the new set, via FW crossover. The asymptotic decay is now 
monotonic. As A is increased further the value of ao of the purely imaginary, n = 0, 
pole then decreases to zero which is equilvalent to the correlation length diverging at 
the spinodal. 



and critical point, with monotonic decay. The phase diagram in figure [T^ shows these 
regions for the mixture with q = 0.5 and A = 0.1. Note that the FW fine again crosses 
the binodal twice, but that these crossings do not occur in coexisting phases. Thus, 
there are coexisting phases which have different types of asymptotic decay, or different 
oscillatory decay wavelengths. 

The most obvious difference that distinguishes the non-additive from the additive 
mixture is the presence of an additional purely imaginary pole in the former case. In 
order to understand where this pole comes from, and to elucidate the full effect of 
introducing non-additivity on the pole structure, we start with an additive mixture and 
slowly increase A from zero. Figure [12] displays the positions of the poles for fixed 
q = 0.5, Tji = ri2 = 0.1 and increasing A from zero. The asymmetric additive hard 
sphere mixture has two sets of complex poles, each set having real components that are 
related to each of the length-scales. As A is increased from zero, a third set of complex 
poles, including the single purely imaginary pole, with very large imaginary components 
appear. As A increases, the imaginary components, Oq, of the new set of poles decrease 
so that this set of poles moves into positions in the complex plane comparable to the set 
of original 'additive' poles. Increasing A further results in the decrease of the imaginary 



Non-additive hard sphere mixtures 



23 




Figure 14. Free interface density profiles, Pi{z), between coexisting fluid phases for 
q = 1 and A = 0.1, and pressures Pa\^/kBT = 3, 5, 7, 9 and 11 corresponding to the 
tie-lines in figure [S] The profiles are plotted as a function of scaled position from the 
interface, z/ an. Since the mixture is symmetric, the density profiles of species 2 are 
identical to those of species 1 under the refiection z — > — z. The insets show the decay 
of \pi{z) — pjl for z < and z > where p\ is the bulk density of species 1 on the 
side of the interface shown in each inset. The profiles in the insets are plotted on a 
semi- logarithmic scale and each profile is offset from the one above by a factor of 10^^. 



parts of the new set of 'non-additive' poles. This proceeds most rapidly for the purely 
imaginary pole, n = 0, (from the new set of 'non-additive' poles) which then becomes the 
leading order pole. For these parameters, the FW crossover occurs at A = 0.101. As A is 
increased further, the purely imaginary pole reaches the real axis, which corresponds to 
the divergence of the correlation wavelength at the spinodal. In the following section, we 
will investigate the repercussion of the asymptotic decay of correlation on the structure 
of the free fluid interface. 



3.3. Structure of the Free Fluid-Fluid Interface 

By numerically minimising fi[pi,p2], using the method outlined in section [23] and the 



planar weight functions given in Appendix C we calculate the equilibrium one-body 
density profiles, Pi(r) (referred to as Pi(r) in the following) for the fluid- fluid interface 
between coexisting phases with a simple iterative Picard scheme. Figure [T3] shows 
the density profiles, Pi{z), of the free interface as a function of distance z from the 
interface for a range of pressures, corresponding to the tie-lines shown in figure |9l Since 
the mixture is symmetric, we plot all results for species 1, but only one representative 
profile for species 2 as an illustration. Starting with coexisting phases close to the critical 
point, we find that the density profiles vary monotonically as a function of z. This is fully 
consistent with the type of asymptotic decay of the pair correlation functions, which is 
monotonic for both coexisting statepoints. As one moves away from the critical point 
one finds that pi{z) becomes oscillatory on the side of the interface where species 1 is 
the majority component {z < 0). This agrees with our results for the asymptotic decay 
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-15 -10 -5 00 5 10 15 

Figure 15. Same as figure [TH but for q = 0.5 and A — 0.1. The profiles for the smaher 
particles are plotted in (a) and those for the larger, species 2, in (b). The coexisting 
pressures are Pafi/kBT — 1.0, 1.25, 1.5, and 1.75 corresponding to the tie-lines in 
figure [121 The insets show the asymptotic decay of the profiles on a semi- logarithmic 
scale, where each profile is oflfset by a factor 10^^ from the one above. 



changing on crossing the FW hne, but there is apparently no oscillations on the side 
of the interface where species 1 is the minority component {z > 0). Similarly, p2{z) is 
oscillatory at Pafi/ksT =11 for z > but does not appear to be oscillatory for z < 0. 
It is clear that as the pressure is increased, the oscillations in the density profiles for 
each species appear on one side of the interface, but to examine how the abrupt change 
in the type of asymptotic decay on crossing the FW line affects the profiles, we must 
investigate the decay of the profiles away from the interface. The inset in figure UM shows 
the intermediate decay of pi{z) for 2; < and for z > 0. We plot the absolute difference 
of the density profiles and their bulk value on that side of the interface, \pi{z) — p^l, on a 
logarithmic scale. The profile between the coexisting phases closest to the critical point, 
Pcxfi/{kBT) = 3, clearly decays monotonically on both sides of the interface. The case 
Pcxfi/{kBT) = 5, appears to be monotonic on the linear plot, but if one looks at the 
intermediate decay behaviour (shown in the inset), one finds that there is oscillatory 
decay on both sides of the interface, but that these oscillations do not appear until at 
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least a distance z/au c^ 3 from the interface. As the pressure of the coexisting phases 
is increased, the oscillations shown in the insets grow in relative amplitude and start to 
appear closer to the interface and thus become more pronounced on the linear scale. 

For q = 0.5 and A = 0.1, corresponding to the phase diagram shown in figure \T2\ 
we again start at the critical point and trace pairs of coexisting state-points along the 
binodal. Figure [15] displays the density profiles for coexisting phases, corresponding to 
the tie-lines in figure [T^ Close to the critical point the coexisting state-points both reside 
in the region of the phase diagram where the asymptotic decay of gij{r) is monotonic. 
On the linear plot we find that the density profiles for coexisting phases close to the 
critical point appear monotonic. The inset of figure [TST a) shows that the density profile 
corresponding to the lowest pressure, Pali/ksT = 1, decays monotonically on both 
sides of the interface. As we increase the pressure and move along the coexistence 
curve, we find that the state-points rich in species 2 crosses the FW line and moves into 
the oscillatory region (labelled Oscillatory-2), while the other state-point remains in the 
monotonic region. Therefore, the density profiles decay with an oscillatory component 
on the side of the interface where species 2 is the majority component {z > 0). The 
inset shows that the profiles for Pafi/ksT = 1.25 and 1.5 both exhibit this behaviour; 
for z < the decay is monotonic and for 2; > the decay is oscillatory. If the pressure is 
increased further, the other state-point (rich in species 1) crosses the FW line and moves 
into the other oscillatory region (Oscillatory-1). The inset shows that the intermediate 
decay of the profiles for Paf^/ksT = 1.75 is oscillatory, but that this is very far away 
from the interface and occurs with a small relative amplitude. These coexisting phases 
both have oscillatory decay but with different wavelengths; if one examines the profile 
for Pa^ /ksT = 1.75, the oscillatory wavelength for 2; > is approximately twice as 
large as that for z < 0. 

3.4- Interface Tension of the fluid-fluid interface 

From the density profiles we have calculated the surface tension of the free interface, 

J=in[p^,p2] + PV)/A, (31) 

where Q[pi{r) , p2{r)] is the grand potential of the inhomogeneous system with the free 
interface, —PV = fi[p\, P2] is the grand potential of the uniform system, and A is the 
area of the interface. Figure [16] displays the surface tension for the mixture with g = 1 
and varying A, plotted against the order parameter \rii — ri^\, where 77^ is the packing 
fraction of species 1 in phase A (and similarly for B). As the mixture is symmetric, 
this quantity is symmetric w.r.t. interchange of species, i.e. \ri^ ~ Vi\ = \V2 ~ V2\- ^^ 
find that increasing the non-additivity has a dramatic effect on the surface tension, for 
constant \ri^ ^ Vi\ = O-l; 7 increases over fifty times between A = 0.1 and A = 1. 

It can be shown [57] that as the critical point is approached, 7 follows a simple 
mean-field scaling law, 7 oc \rii — rj^l^. In order to check our calculations of 7 we plot 7 
against \r]f^ — ri^\ on a double logarithmic scale in the inset of figure [THl For comparison, 
in the inset we show the asymptotic result, 7 = a\r]^ — r]f\^ (labelled x^) where a is 
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Figure 16. The surface tension, 7, of the planar fluid-fluid interfaces with parameters 
q = 1 and A — 0.1, 0.2, 0.3, 0.4, 0.5 and 1, plotted against the absolute difference in 
the partial packing fraction of species 1, |?yi — ?7f |, in the two coexisting phases A and 
B. The inset shows the same quantities on a double logarithmic scale, and compares 
them to the mean-field behaviour, 7 oc |?7^ — ?7f P, labelled x^ (thin solid line). 
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Figure 17. Same as figure [TBI but for q = 0.5 and A = 0.1, 0.2, 0.4 and 0.5, and 
shown as a function of the absolute difference in the total packing fraction, [77 — jy^l, 
in the two coexisting phases A and B. Within this representation we find that for fixed 
order parameter, |?7 — Jy^l, 7 varies non-monotonically with A. The inset shows the 
same quantities on a double logarithmic scale and compares them to the mean-field 
behaviour, 7 oc \ri^ — rp\^ ^ again labelled x^ . 



a proportionality constant. For all values of A, as \ri^ — rjf\ approaches zero, 7 tends 
towards the mean-field behaviour, i.e. the slope of the curves in the inset tends towards 
the slope of the asymptotic result. 

Figure [17] displays the surface tension of the free interface for the mixture with 
q = 0.5 and varying A. As the mixture is asymmetric, we plot these results using the 
(species independent) order parameter \ri^ — ri^\ which is the absolute difference in the 
packing fraction between phases A and B. Note that 7 is scaled by the square of the 
diameter of the larger species, (J22, so as to keep the values of 7 comparable across a 
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Figure 18. Same as figure [T71 but for q = 0.1 and A =0.2, 0.4 and 1. 

range of q values. In this representation, we find that the curves are broadly similar 
to those in figure Uni but that 7 exhibits a rapid increase with increasing \rj^ — r]^\. 
Furthermore, for all values of the order parameter, the value of 7 for A = 0.2 is smaller 
than the value for A = 0.1. These two features arise from our choice of order parameter. 
If we use an order parameter similar to the one in figure [161 we do not have either of these 
two features. Although we have a different order parameter, the surface tension follows 
a similar scaling law as we approach the critical point, 7 oc \rj^ — rj^\^. In figure ITS] we 



plot 7 for the mixture with fixed g = 0.1 and varying A as a function of |?7"^ ~ V^l- We 
find that, unlike in figure [T71 7 does not increase rapidly as |?7"^ — ?7^| approaches its 
maximum value, and that 7 increases monotonically with A (for the values considered). 



4. Discussion 



Using a fundamental measure density functional theory we have investigated some of the 
properties of homogeneous and inhomogeneous fiuid states of a binary non-additive hard 
sphere model with positive non-additivity. This model exhibits fiuid-fiuid demixing. We 
have calculated the coexistence curves and showed that these compare reasonably well 
to existing simulation results. The theory predicts that the critical point occurs at a 
pressure and density lower than the simulation results. This is typical of mean-field type 
theories, such as DFT, which do not take account of all fiuctuations in the fiuid. We 
have not investigated whether the fiuid-fiuid phase transitions are stable with respect 
to crystallisation, which is expected to occur at high packing fractions. To investigate 
this one would require a more sophisticated functional which is capable of modelling 
the extreme confinement in a crystal. Moreover, even for the additive mixture, where a 
suitable theory exists [SHI EO] , we are not aware of any systematic DFT investigation of 
freezing. 

Using the Ornstein-Zernike equation, we calculated the asymptotic decay of 
correlation functions, Qijir), by solving for the poles of the partial structure factors. 
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Sij{k), in the complex plane. Using Cauchy's theorem, one can express the correlation 
functions as an infinite sum over these complex poles. In particular the poles with the 
smallest imaginary part are interesting, as these determine the asymptotic, r — > oo, 
decay of the entire set of gij{r). We find that for q < 1 and A > there is, at low 
densities, a crossover between two modes of asymptotic damped oscillatory decay with 
different wavelengths, which are similar to the diameters of the two species. For A > 
we find Fisher- Widom crossover from oscillatory to monotonic asymptotic decay as the 
coexistence region is approached. We find that the positive non-additivity introduces a 
new set of poles, including one purely imaginary pole. As A is increased from zero this 
new set of complex poles appear in the complex plane initially with very large imaginary 
components. As A is increased, the value of the imaginary components decreases and 
the poles occupy a region of the complex plane similar to the set of poles that exist 
already in the additive model. 

One might imagine that the new length-scale, ai2, would induce a third regime 
where the asymptotic decay is oscillatory with a wavelength similar to the cross-species 
diameter, ai2. However, for A > this does not occur since the new set of complex 
poles, related to a non-zero Ru, is always accompanied by a purely imaginary pole 
which is always the leading order pole of this set. It would be interesting to investigate 
the case A < in future work. 

Furthermore, we have studied the inhomogeneous free fiuid interface between 
coexisting phases and have calculated the density profiles and the surface tension. We 
showed how the type of asymptotic decay affects the intermediate and short-range 
behaviour of the density profiles. We have presented detailed results for the surface 
tension of the free fiuid interface. These can be compared to both simulation and 
experimental results, and furthermore play a vital role in the investigation of capillary 
condensation phenomena. 

Appendix A. Weight Functions in Fourier Space 

For completeness we include the Fourier space representations of the weight functions. 





w^ = A7r{s-kRc)/k^, 


W2 


= AnRs/k, 




wi = {kRc + s)/{2k), 


Wo 


= c+{kRs/2), 


and 










w\ = 4:7r{kRc + s)/k, 




wl = c-{kRs/2), 




wU = -^{k'Rc + 3ks), 




w_i = {k^Rc-ks)/2, 




W-2 = —-r^k^Rs, 




w.s = ^(k'Rc-Ses), 


where s = 


= sin(A;i?) and c = cos{kR). 







Appendix B. Free Energy Density Contributions 

Representative cases of the free energy terms , ^a/s, where, a, f3 = to 3, are shown 
below. The remaining terms can be obtained through symmetry by changing the species 
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labels: ^pa = ^ap{nv ^ ni ,nr ^ Ur )■ rj is the the total packing fraction, given by 



(1) , (2) 



3 / ,inN 3 ,^. ,^, ,-,. / ,^.N 3 



nf ^n« (n?)) (n«) nf^nf 41)^(2)^2) ^(1)^(1)^(2) n« (n^^) 

•^00 — -4 1 -4 1 ^2 ' ' ^^5 1- 



47r(l-r/)* 47i{l-riy (1-^) (1-^) 127r (1 - r;)' 

^ (1-^)^ ^ 127r (1 - r^)^ ^ 2An^l-r]f ^ 1 - r/ ' 



+ 



\4 ' 



327r2 (1 - T]) 

247r(l-r7) 1-^ 

nr'nr "1 V"'2 ; ("'2 ) ^1 (^2 ) ("-2 



" 1-^ ^ Stt (1 - 77)' ^ Stt (1 - r^)' ^ 327r2(l-r/)' ' 



(1) (2) 

$22 = ' 

\ — rj 

•^23 = -4^^1n(l-r7), 

$33 = (l-?7)ln(l-r7)+r]. 

Appendix C. Weight Functions in Planar Geometry 

In this paper we consider planar density profiles, where we can simplify the convolutions 
by performing the integration over the radial direction in advance, yielding a set of planar 
weight functions; 



w^^Hz) 



27ry""da4^^(\/e+^), (c.i) 



where ^ = y a;^ + y"^. This yields 

m{z) = 7rsgn{R)e{R - \z\){\R\^ - \z\^), 
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and 



W2{z) = 2nRQ{\R\ - \z\), 

w,{z) = ^ sgn(i?)[0(|i?| - 1^1) + z5i\R\ - \z\)], 

Wo{z) = ^-6{\R\-\z\)-y6'{\R\-\z\), 



wliz) = 27isgn{R)[Q{\R\ - \z\) + z6{\R\ - \z\)], 



Wo[Z) 



\R\-\z\) + z6'{\R\-\z\ 



wU^) = :^[6'i\R\-\z\) + z5^'\\R\-\z\ 



w^i[z) 
w-2{z) 

W^3{z) 



-sgniR)5i\R\-\z\) + z5^'\\R\-\z\)], 
[36^^\\R\ - \z\) - z6'^'^\\R\ - \z\)], 



327r 



±-sgn{R)[-76^'\\R\ - \z\) + z6^'H\R\ - 



Furthermore, we can perform the convolution of a general one-dimensional function, 
f{z), with these planar weight functions. 



F^^\z)= dz'f{z')w^^\\z-z'\), 



givmg 



pz' + \R\ 

Fs{z) =nsgn{R) / dz'f{z')[R' 

Jz'-\R\ 

F2{z) = 27rR / dz'fiz'), 

Jz'-~\R\ 

"z' + lR 



[Z — Z 
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oZn 



. ± 



±/'="(:±l-R|) 



+ 5^E/'^'(^±W) 



327r 

where f'{z), f"{z), /'■'^•'(-z) and f^^\z) represent successive derivatives of f{z). 

The weighted densities and convolutions with weight functions for the individual 
species are calculated using Fourier transforms, while the convolutions with the kernel 
functions are calculated directly by integrating over z'. The 7th derivative of the 
free energy density derivative term, d$a/3/d^7, is calculated using a central difference 
approximation with a symmetric (7 + l)-point stencil. 
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